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PREFACE 
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Mr.  D.  E.  Howell,  Chief,  Information  Sciences  Division,  and  Mr.  L.  A.  Gambino, 
Director,  Computer  Sciences  Laboratory. 

A  portion  of  the  computer  implementation  was  performed  by  Cadet  Glenn 
All ton  of  the  Air  Force  Academy  and  by  Mr.  William  Moore,  a  student  at 
Virginia  Polytechnic  Institute. 

COL  Edward  K.  Wintz,  CE,  was  Commander  and  Director  and  Mr.  Walter  E.  Boge 
was  Technical  Director  of  the  Engineer  Topographic  Laboratories  during  the 
study  preparation. 
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REGISTRATION  OF  A  LANDSAT  IMAGE  TO  A 
DIGITAL  TERRAIN  MATRIX— AN  ERROR  ANALYSIS 


INTRODUCTION 


^One  reason  for  registering  a  Landsat  digital  image  to  a  Digital' 
Tterr»i»- Matrix  {DTM)4>is  to  add  another  dimension  to  the  signature  used  in 
scene  classification.  Ground  elevation  associated  with  a  pixel  can  be 
converted  to  a  terrain  slope  that,  along  with  orientation  with  respect  to 
the  sun,  can  be  used  as  an  additional  component  in  a  pattern  recognition 
exercise.  Pixel  heights  can  also  be  used  to  restrict  the  pattern 
recognition  scheme  to  test  over  the  most  likely  subset  of  classes  when 
scene  classes  can  be  organized  by  terrain  heights  a  priori.  The  purpose 
of  this  report  is  to  evaluate  the  accuracy  with  which  the  elevation 
assignments  can  be  made  when  the  registration  process  described  herein  is 
used. ^ 

RELEVANT  GEOMETRY 

The  dynamic  Landsat  imaging  event  will  be  modeled  by  requiring  that 
the  exterior  orientation  parameters  of  the  camera  system  be  adjusted  to 
fit  common  data  points  identified  on  the  DTM  and  on  the  Landsat  digital 
image.  These  parameters  include  the  time-dependent  vehicle  position  and 
the  time-dependent  vehicle  coordinate  frame  attitude.  The  Landsat 
interior  orientation  parameters  will  not  be  adjusted.  These  parameters 
include  the  detector  sizes,  the  camera  focal  length,  scan  angle  range, 
scan  angle  cycle  time,  and  scan  time.  The  geometry  that  relates  Landsat 
coordinates  (e  =  line  and  s  =  sample)  to  ground  coordinates  defined  by  the 
DTM  is  given  in  detail  in  appendix  A. 

Procedure.  There  are  at  least  two  ways  to  identify  common  points  in 
the  DTM  and  in  the  Landsat  digital  image.  In  the  first  case,  common 
points  are  identified  on  maps  and  on  Landsat.  The  common  data  points  will 
be  used  to  adjust  the  Landsat  geometric  model.  If  the  point  identifica¬ 
tion  is  accurate  and  if  there  is  negligible  bias  between  the  DTM  and  the 
map,  then  the  intersection  procedure  described  in  appendix  A  will  produce 
accurate  results. 

The  second  procedure  relates  the  DTM  and  the  Landsat  digital  image 
directly  by  viewing  both  at  the  same  time  on  separate  TV  displays.  A  gray 
shade  relief  image  is  used  to  transform  the  DTM  into  a  terrain  surface 
that  is  recognizable  as  such  (see  appendix  E).  The  analyst  indicates 
common  points  on  the  two  TV  displays  with  cursors.  Note  that  in  either 
pointing,  the  analyst  is  in  direct  communciation  with  the  pertinent  data 
base. 


-  V- *.>*.*■ s, ■.  .  •  *vVv  ■  /vv-.-v 

'v’v'Vv*- \ 
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Adjustment  Model.  At  any  Instant  of  time,  there  are  six  exterior 
orientation  parameters;  namely,  the  three  angles  that  define  the  attitude 
matrix  of  the  DTM  reference  frame  with  respect  to  the  vehicle  frame  (Mqv  ) 
and  the  three  coordinates  that  locate  the  vehicle  origin  in  the  DTM.  The 
three  fundamental  angles  that  define  will  be  modeled  as  time 
polynomials.  For  example,  if  the  fundamental  angles  are  roll  (w) ,  pitch 
($  ),  and  yaw  (K),  then 

©  t  ‘  B  T» 

where 


B:  3X(N+1)  matrix  of  attitude  constants 

V =  o»  t*  t2»  •••»  tN) 


Vehicle  positions  will  also  be  modeled  by  time  polynomials. 


where 


A:  3X(MH)  matrix  of  position  constants 
TM  =  (1»  t,  t2,  ...»  tM) 


The  following  discussion  describes  how  second-order  time  polynomials 
provide  sufficient  accuracy  for  each  of  the  six  exterior  parameters.  In 
fact,  the  three  attitude  parameters  can  be  represented  by  linear  time 
polynomials  with  no  more  than  one-half  sample  spacing  error. 

The  actual  Landsat  orbit  was  assumed  to  follow  a  Keplerian  path,  and  a 
30-second  ephemeris  in  earth-fixed  coordinates  of  such  an  orbit  was 
developed  using  the  parameter  values  in  table  1. 


Table  1.  Orbital  parameters 

Tq  «  0.0  :  Time  of  Perigee  Passage 

o 

=  225  :  Longitude  of  Ascending  Node 

o 

1  =  80.79°  :  Inclination 

w  =  135°  :  Argument  of  Perigee 

e  =  0.001019  :  Eccentricity 

a  =  7285989  meters  :  Semi-Major  Axis 


Note  that  ^NQ  is  measured  westward  from  Greenwich  and  pertains  to 
time  T0.  The  inclination  is  measured  counterclockwise  at  the  ascending 
node.  The  earth  is  regarded  as  a  sphere  with  radius,  R  =  6,359,550 
meters,  which  means  that  the  Landsat  vehicle  height  is  H  £  (a-R)  =  926,439 
meters.  A  30-second  pass  corresponds  to  a  193.7-kilometer,  along-track 
image,  which  is  larger  than  the  conventional  picture. 

The  vehicle  reference  frame  is  nominally  oriented  so  that  Zy  is  normal 
to  the  surface  of  the  earth  and  so  that  Xy  is  tangent  to  the  flight 
path.  The  Yy-axis  completes  a  right-handed  system.  A  set  of  roll,  pitch, 
and  yaw  angles  was  generated  for  each  of  the  31  ephemeris  positions. 

These  data  as  well  as  the  position  data  were  fit  in  turn  to  linear, 
quadratic,  cubic,  and  finally  to  quartic  time  polynomials.  The  standard 
deviations  of  the  residuals  were  calculated  for  each  of  the  six  exterior 
orientation  values  and  for  each  of  the  four  polynominal  models.  These 
values  represent  discrepancies  between  the  Keplerian  model  and  the 
polynomial  models  for  a  30-second  Landsat  exposure.  The  standard 
deviations  are  given  in  table  2.  The  position  discrepancies  are  expressed 
in  meters  and  the  attitude  discrepancies  are  expressed  in  seconds  of  arc. 

The  total  attitude  error  at  Landsat  heights  (less  than  106  meters)  is 
less  than  0.1  meter  per  point  for  the  qudratic  model,  and  the  total 
position  error  is  less  than  0.7  meter.  The  corresponding  errors  for  the 
linear  model  are  36  meters  and  274  meters,  respectively.  The  linear  model 
error  for  position  is  large  when  compared  to  the  Landsat  foot  print 
dimension  of  79  meters. 
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Table  2.  Model  discrepancies 


Exterior  Orientation 


Model 

X 

Y 

Z 

V 

♦ 

K 

Linear 

1.4X102 

1.4X102 

1.9X102 

4.39X10° 

2.98X10° 

5.17X10° 

Quadratic 

3.4X10-1 

3.8X1CT1 

5.3X10"1 

1.4X10*2 

Osi 

1 

C 

r—i 

X 

CO 

• 

r—( 

8.2X10-3 

Cubic 

8.4X10"4 

6.6X10-4 

1.0X10-3 

3.1X10"5 

3.0X10-5 

1.9X10-5 

Quartic 

4.1X10-4 

1.6X10"4 

4.1X10'4 

9.7X10-6 

4.5X10-6 

5.3X10'6 

ERROR  ANALYSIS 

The  30-second  Landsat  exposure  described  above  along  with  the  data 
generator  described  in  appendix  C  was  used  to  develop  28  sets  of  exact 
data.  For  convenience,  the  points  were  defined  in  a  local  frame  centered 
at  t  *  15  seconds.  The  points  are  displayed  in  figure  1  along  with  the 
Landsat  ground  track;  the  grid  spacing  in  figure  1  is  20  kilometers.  Note 
that  the  ground  tract  is  from  west  to  east.  The  reason  for  this  is  that 
the  inclination  should  have  been  (ir-i)  =  99°21,  instead  of  i  =  80?79. 
Except  for  earth  rotation,  the  two  situations  are  exactly  symmetric  about 
the  pole,  which  means  that  using  the  supplement  of  the  inclination  will 
have  no  effect  on  the  error  analysis  of  this  report. 

Least-Squares  Adjustment.  Five  subsets  of  the  28  ground  points  were 
selected  to  represent  five  well-distributed  control  point  patterns  for 
adjustment.  In  addition,  four  sets  of  exposure  station  weights  were 
processed  to  determine  the  effect  of  a  priori  orbital  data  on  the 
adjustment.  The  a  priori  exposure  station  standard  deviations  are  given 
in  table  3. 


Table  3.  Land sat  epheneric  weights 

Height 

qc(m)  qc(m/sec)  °c(m/sec2) 


W1 

100 

1.0 

0.1 

W2 

500 

2.0 

0.2 

W3 

1000 

5.0 

0.5 

Wa 

2000 

10.0 

1.0 

The  five  control  point  patterns  are  defined  in  table  4  and  are 
referenced  to  figure  1. 


Table  4. 

Control  point  patterns 

Control  Point 
Pattern 

Number 

Points 

C1 

28 

All  28 

C2 

11 

#1,  #2,  #4,  #7,  #8,  #12,  #13,  #16 
#19,  #24,  #28 

c3 

6 

#1,  #2,  #4,  #16,  #19,  #28 

C4 

5 

#1,  #4,  #16,  #19,  #24 

C5 

17 

#1,  #3,  #4,  #5,  #6,  #7,  #8,  #9, 
#10,  #11,  #12,  #13,  #14,  #15,  #16 
#17,  #19 

The  data  in  tables  3  and  4  were  used  to  develop  the  several  18  by  18 
covariance  matrices  of  the  exterior  orientation  parameter  estimates. 


Table  5.  Standard  errors  In  pixels  for  Wj 
Points 


Pattern 

C1 

c2 

C3 

C4 

C5 


2 

2 

2 

4_ 

2 

3.4 

3.4 

4.1 

2.9 

2.1 

0.8 

0.8 

0.9 

0.5 

0.6 

5.0 

5.0 

5.9 

3.7 

3.3 

1.0 

1.0 

1.2 

0.6 

0.7 

12.6 

12.5 

10.6 

6.5 

31.8 

1.3 

1.3 

6.1 

0.6 

1.5 

35.8 

35.6 

151.9 

6.7 

11.5 

10.2 

10.0 

3.3 

1.0 

0.9 

4.1 

4.1 

8.8 

3.9 

2.5 

0.8 

0.8 

1.3 

0.6 

0.6 

2 

1-2 

1-3 

4-5 

4-6 

3.3 

0.01 

5.3 

3.9 

4.0 

0.8 

0.01 

1.0 

0.6 

0.7 

4.8 

0.01 

7.6 

5.4 

5.8 

0.9 

0.01 

1.4 

0.7 

0.9 

6.5 

0.02 

21.2 

37.3 

10.6 

0.9 

0.02 

5.6 

1.8 

1.0 

178.8 

0.16 

116.8 

10.9 

174.7 

3.4 

0.14 

7.6 

1.1 

2.8 

6.9 

0.01 

9.9 

5.3 

6.8 

0.9 

0.01 

1.4 

0.7 

0.9 

Table  6.  Standard  errors  in  pixels  for  W2 
Points 


Pattern 


C1 

c2 

c3 

C4 

C5 


2 

2 

2 

2 

2 

3.5 

3.4 

4.2 

3.2 

2.1 

1.1 

1.1 

1.8 

1.0 

1.4 

5.1 

5.0 

6.1 

4.1 

3.4 

1.3 

1.3 

2.2 

1.1 

1.7 

13.8 

13.8 

11.4 

6.5 

35.6 

1.6 

1.6 

6.5 

1.4 

1.8 

41.6 

41.4 

176.2 

6.8 

13.2 

11.6 

11.4 

3.9 

1.5 

1.9 

4.1 

4.1 

9.0 

4.5 

2.6 

1.1 

1.1 

2.3 

1.1 

1.7 

2 

1-2 

1-3 

4-5 

4-6 

3.4 

0.01 

5.6 

4.3 

4.3 

1.7 

0.01 

1.3 

0.7 

1.1 

4.9 

0.01 

7.9 

5.7 

6.2 

2.0 

0.01 

1.8 

1.0 

1.4 

6.6 

0.02 

23.5 

40.8 

10.7 

2.1 

0.02 

6.2 

1.8 

1.4 

208.1 

0.19 

135.2 

12.6 

204.4 

4.6 

0.16 

8.8 

1.4 

3.5 

7.0 

0.01 

10.0 

5.9 

7.2 

2.0 

0.01 

1.8 

1.0 

1.4 

M 
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DISCUSSION 


The  expected  DTM  horizontal  standard  errors  owing  to  imprecision  in 
the  exterior  orientation  estimates  are  summarized  in  figures  2  through 
6.  Note  that  the  results  are  presented  in  terms  of  the  Landsat  foot  print 
dimension  (79  meters)  and  that  the  results  pertain  to  a  standard  error  of 
one  line  and  one  sample  in  the  DTM  point  identification. 

Except  for  the  results,  the  along-track  errors  (cre)  do  not  vary 
greatly  as  the  exterior  orientation  position  parameter  weights  vary.  The 
across-track  errors  (og)  do  vary,  but  are  small  with  respect  to  the  along- 
track  errors.  One  reason  for  the  noticeable  increase  in  error  for  as 
the  weights  decrease  is  that  there  was  only  1  degree  of  freedom  associated 
with  the  C,  adjustment.  Patterns  Cp  Cj,  and  are  similar  in  that  the 
control  points  are  symmetrically  distributed  over  the  Landsat  format.  The 
along-track  errors  for  Cj  and  C2  approximate  the  /""N  law  of  least-square 
averages;  whereas,  there  is  little  difference  in  the  across-track  error. 
The  error  results  increase  considerably  from  C2  because  there  are  only 
3  degrees  of  freedom  associated  with  the  adjustment.  Finally,  there  are 
more  degrees  of  freedom  associated  with  (25)  than  with  C2  (13),  yet  the 
error  estimates  are  larger.  The  reason  for  this  is  that  there  were  no 
observations  along  the  vehicle  ground  track  for  C^. 

In  order  to  lower  the  horizontal  intersection  errors,  either  more 
points  must  be  identified  and  measured  or  independent  observations  on  the 
attitude  parameters  must  be  obtained.  A  considerable  effort  must  be 
expended  if  the  horizontal  errors  are  to  be  reduced  without  some 
additional  control  in  the  attitude  parameters.  Let  a  be  the  measuring 
and  Identification  Landsat  error,  expressed  In  pixels™  associated  with  the 
initial  DTM  point  identification.  Using  the  Cj  results  as  a  base,  the 
following  approximation  can  be  used  to  estimate  the  expected  standard 
error  in  the  along-track  horizontal  intersection  error: 


°M  <  °e  < 


N  =»  Number  of  control  points. 


For  example,  if  N=95  and  if  the  points  are  symmetrically  scattered 
over  the  Landsat  format,  then  om  <  ae  <  2am,  i.e.,  the  along-track 
intersection  standard  error  is  between  80  and  160  meters.  This  result 
pertains  to  a  point  identification  and  measuring  error  of  1  pixel,  which 
is  unrealistic.  Note  that  identification  errors  can  easily  assign  the 
measurement  to  the  next  swath,  which  is  a  time  error  of  0.07342  second  or 
6  pixels. 
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Figure  4.  Error  suwary  for  C, 
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Figure  5.  Error  suaaary  for  C^. 
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Figure  6.  Error  suanaty  for  C*. 


The  technique  described  above  and  in  appendix  E  is  one  way  to  relate  a 
Landsat  image  to  a  DTM.  It  is  difficult  to  believe  that  enough  points  can 
be  accurately  defined  and  measured  on  both  of  the  soft  copy  representa¬ 
tions  to  precisely  determine  the  exterior  orientation  of  the  scanning 
vehicle  in  the  DTM  coordinate  frame.  Consider  figures  7  and  8.  Figure  7 
represents  a  gray  shade  subscene  derived  from  a  DTM,  which  in  turn  was 
produced  by  DMA  from  1:250,000  scale  maps.  The  horizontal  ground  spacing 
of  the  DTM  is  64  meters.  Figure  8  is  a  corresponding  Landsat  subscene. 

It  seems  unlikely  that  the  procedure  will  work  in  this  situation.  It  is 
true  that  the  area  is  relatively  flat.  Figure  9  is  another  gray  shade 
relief  (and  of  another  region)  where  topographic  detail  is  more  obvious 
and  more  readily  identified. 

If  points  are  identified  on  maps  and  on  Landsat  and  if  a  DTM  is 
intersected  to  complete  the  rectification,  then  there  will  be  the 
additional  concern  of  the  distortion  and  imprecision  between  the  DTM  and 
the  map.  Considering  the  expense  and  likelihood  of  an  accurate 
registration,  perhaps  it  should  be  shown  beforehand  that  the  addition  of 
terrain  elevation  information  to  the  pattern  recognition  exercise  will  be 
worth  the  effort. 


It  appears  that  the  full  point-by-point  rectification  as  described  in 
appendix  A  will  be  an  expensive  and  time-consuming  effort,  especially  if 
the  absolute  and  relative  errors  are  to  be  held  to  about  5  pixels  (400 
meters).  Perhaps  a  more  practical  approach  would  be  to  intersect  the  DTM 
using  a  regular  pattern  of  Landsat  pixels  and  then  to  contour  the 
resultant  pattern,  which  could  then  be  easily  displayed  on  the  Landsat 
soft  copy. 


CONCLUSIONS 

1.  The  exterior  orientation  of  a  Landsat  exposure  can  be  accurately 
represented  by  quadratic  time  polynomials  for  each  of  the  six  parameters. 

2.  The  technique  of  determining  corresponding  points  between  a 
Landsat  image  and  a  gray  shade  relief  image  is  not  feasible  over  the 
relatively  flat  terrain  used  in  this  study. 

3.  The  utility  of  terrain  elevation  information  in  pattern 
recognition  should  be  demonstrated  before  the  rectification  procedure 
described  in  this  report  is  used. 
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1=1,6 


Specific  line  within  Kth  swath 


J 


Sample  value  along  a  line 


K 

dITK 


Specific  swath  in  which  P  was  imaged 

Distance  between  P  and  V  at  exposure  time 

Direction  cosines  of  imaging  ray  in  vehicle  reference 
frame  (V) 


M, 


GV 


JK 


Rotation  matrix  relating  (G)  to  (V)  at  exposure  time 


The  swath  notation  is  necessary  since  Landsat  collects  data  simultaneously 
along  six  parallel  lines.  The  subscripts  I,  J,  and  K  define  which  of  the 
basic  measurements  have  a  bearing  on  the  particular  calculation. 

The  camera  for  any  one  of  the  four  multispectral  bands  can  be  repre¬ 
sented  by  a  simple  linear  array  (see  figure  A2)  in  the  focal  plane.  An 
oscillating  flat  mirror  scans  in  such  a  manner  that  imagery  moves  across 
the  array  in  the  y-direction,  while  the  vehicle  carries  the  scanner  along 
in  the  x-direction.  A  line  of  imagery  moves  across  the  four  arrays  where 
a  given  wavelength  is  sampled  at  each  array.  The  detectors  are  sampled  so 
that  the  six  corresponding  points  from  band  to  band  are  equally  spaced  in 
time.  Except  for  a  negligible  image  motion,  the  four  spectral  values  for 
any  one  of  the  six  pixels  pertain  to  the  same  ground  foot  print. 

zs 


The  intersection  model  expressed  by  equation  1  pertains  to  any  one  of 
the  four  spectral  bands  and,  in  fact,  to  a  composite  image  derived  from 
two  or  more  of  the  spectral  bands.  This  is  true  because  the  corresponding 
foot  prints  are,  except  for  image  motion,  identical.  Let  figure  A2 
represent  the  instantaneous  array  of  detectors  for  a  specific  band,  or  for 
a  particular  composite.  The  scanner  coordinates  for  the  I1*1  point  are 


sin  Y 


v-cos  Y 


where 


2  2 
Xj.  +  f 


The  six  discrete  values  can  be  computed  directly  from  the  focal  length 
and  from  the  detector  dimensions.  For  example,  the  focal  length  for 
Landsat  2  is  32.34  inches  and  for  each  detector  is  71  pm  by  71  urn. 

The  remaining  model  parameters  are  described  in  figure  A3.  The 
vehicle  reference  frame  (V)  is  nominally  oriented  so  that  Zv  is  normal  to 
the  earth's  surface  and  so  that  Xy  is  tangent  to  the  flight  path.  The 
three  fundamental  angles  used  to  define  Mgy  will  be  modeled  as  time 
polynomials.  For  example,  if  the  fundamental  angles  are  roll  (w),  pitch 
(<fr),  and  yaw  (K),  then 


where 


tJK*  t2JK»  *'•»  fcNJK^ 

TJK  -  To 
Exposure  time 

Arbitrary  epoch  time 

3X(N+1)  matrix  of  attitude  constants 

Order  of  time  polynomial 


Instantaneous 


Vehicle  positions  will  also  be  represented  by  time  polynomials. 
(*v)  *  ATh 

\  Zv/  JK 
where 

A  :  3X(M+1)  matrix  of  position  constants 

M  :  Order  of  time  polynomial 
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The  direction  cosines  of  the  imaging  ray  in  (V)  are  derived  from  the 
scanner  coordinates  and  from  MSVj  ,  the  rotation  matrix  relating  the 
oscillating  flat  mirror  to  (V).  £rora  figure  A3, 


0 

Cos  £ 
- 


-si  c\ 

Cos  KjJ 


(A5) 


where  £^  is  the  scan  angle. 

Then,  the  direction  cosines  of  the  imaging  array  in  (V)  are 


Cos 

Cos 


(A6) 


The  coordinates  of  P  in  (G)  can  be  calculated  from  a  Landsat  image, 
providing  the  following  quantities  are  known  or  can  be  computed: 


A  and  B 


d 

IJK 


3(M+N+-2)  exterior  orientation  constants 

Exposure  time 

Scan  angle 

Array  measurement 

Slant  range  distance  between  the  vehicle  and  P 
at  exposure  time. 


The  exterior  orientation  constants  are  determined  by  a  least-squares 
adjustment  (see  appendix  B).  Exposure  time  is  computed  as  follows: 


T 


JK 


Tl  +  (K-l)AT  + 


6T 
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:  Time  of  start  of  first  swath 


AT 

K 

AT 

J 


Scan  time,  plus  return  time 
Swath  number 
Scan  time 
Sample  number 

Total  number  of  samples  per  line 


The  scan  angle  is  computed  in  the  following  way: 


where 

£  :  Scan  angle  for  J“1 

r 

?  :  Scan  angle  for  J=NP 

L 

The  following  values  were  used  in  Landsat  2: 

AT  =  0.073A2  seconds 

AT  =  0.033  seconds 

Cp  =  -5.8  degrees 

=  +5.8  degrees 


The  slant  range  distance  d,™  can  be  calculated  only  if  additonal 
information  about  the  terrain  is  provided,  a  DIM  for  example.  Rewrite 
equation  (Al)  and  drop  the  subscripts  I,  J,  and  K, 


\'‘p>  ~p>  Up/ 


P*  P*  'P'  GV„ 


Suppose  a  good  estimate  of  Zp,  say  Zp  ,  exists. 


then  from  equation  (A7) 


and 


(Zp^  -  Zy)up 


=  Xu  +  d_A 


o  P 


=  Y„  +  d  tn 
V  op 


An  improved  value  of  Zp  can  be  obtained  by  interpolating  in  the  DTM. 
Consider  figure  A4. 


Figure  A4.  Bilinear  interpolation. 


Bilinear  interpolation  for  zPj  produces  the  following  values: 

zPj  =  ZA  +  6x(ZB  -  Zh> 

ZA  =  Z(i,j)  +  6Y[Z(i,j+l)  -  Z( i , j)  ] 

ZB  -  Z( i+1 , j )  +  5y[Z( i+1 , j+1)  -  Z(i+l,j)] 

6X  =  AX/ A 

6Y  =  AY/ A 

A  :  Horizontal  spacing  of  DTM 

If  |ZPQ  -  zPi|i  1,  a  preselected  value,  then  Zp  =  zPj.  If  the  convergence 
criterion  is  not  met,  then  repeat  the  process  using  dj  =  (zpj  -  Zy)/ u  and 
associated  values  for  ^Pj  and  ^p^.  ^ 

Line  of  sight  should  not  be  a  problem  here.  Consider  figure  A5. 


Z 


Figure  A5.  Line  of  sight. 


Either  Pj,  ?2 ,  or  Pg  would  satisfy  the  iterative  process,  depending  on  the 
^p|  initial  estimate  Zp  .  From  the  diagram  Tan  £  =  AY/AZ,  or 


AZ  =  AY/Tan  £  >  9.845  AY 


Since  |  £  |  <  5.8  degrees 


Line  of  sight  is  obscured  if  the  change  in  elevation  AZ  exceeds  9.845  Ay. 
The  minimum  Y-spacing  for  a  DTM  derived  from  1:250,000  scale  maps  is  208.3 
feet  (65  meters),  which  says  that  AZ  must  change  by  at  least  2051  feet 
(650  meters)  over  the  horizontal  spacing  before  line  of  sight  fails. 

The  number  of  iterations  should  be  monitored  since  it  is  possible,  but 
unlikely,  that  the  process  could  cycle.  Consider  figure  A6. 

Z 


Figure  A6.  Iterative  breakdown. 


The  process  will  cycle  if  the  slope  at  P  and  along  1  is  S  =  90°  -  £.  The 
linear  section  1  is  defined  by  the  DTM  and  by  the  plane  containing  the 
local  vertical  at  V  and  along  the  ray  D.  Since  |  £  |  <  5.8  degrees,  S  must 
be  >  84.2  degrees  for  the  process  to  cycle. 


31 


.V_\  .V.VA 


APPENDIX  B.  LEAST- SQUARES  ADJUSTMENT  OF  LAND SAT 
EXTERIOR  ORIENTATION  MODEL  PARAMETERS 


i 


Consider  equation  (Al).  Drop  the  subscripts  I,  J,  and  K,  and  rewrite 
the  equation  as  follows: 


M 


xp-  xv 


GVr 


Y  -  Y 
P  V 


1/d 


Z  -  Z 
P  v. 


PV 


Perform  the  indicated  matrix  operations,  and  divide  the  resulting  first 
and  second  equations  by  the  third  equation  to  obtain  the  following 
condition  equations: 


F  =  Tan  Y  Sec  5  + 


“11<VV  +  m12<W  +  m13(VV 


=  0 


(Bl) 


m31<Xn-X„>  +  m32<Yn  -Y,)  +  m33^Xr.~^t») 


P  V 


P  V 


P  V' 


G  *  Tan  £ 


+  n>2l(Xp-Xv)  +  m22(Yp-Yv)  +  m23(Zp-Zv) 


m3l(XD~X„)  +  ®32(Yn-Y„)  +  «33(Z«-Z„) 


p  -V'  J32^p-V  33'"P  "V' 


0 

0 


(B2) 


where  (m^j;  i=l»3  and  j=l,3)  are  the  elements  of  MqV 


The  expressions  F  and  G  are  nonlinear  in  the  exterior  orientation 
constants.  In  order  to  use  the  least-squares  technique,  F  and  G  must  be 
approximated  with  linear  equations.  Assume  that  reasonable  estimates  of 
the  3(N+FH-2)  unknown  constants  defined  in  equations  (A3)  and  (A4)  exist. 
Equations  (Bl)  and  (B2)  will  be  approximated  by  the  constant  and  first- 
order  Taylor  expansion  terms  of  F  and  G  about  the  estimates.  The 
approximate  equations  are  linear  in  corrections  to  the  parameter 
estimates,  and  thus  the  least-squares  technique  can  be  used  to  solve  for 
the  corrections. 
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Represent  the  condition  equations  (Bl)  and  (B2)  from  equation  (Al)  as 

F  =  m  +  R/T  =  0 
G  =  v  +  S/T  =  0 


A  straightforward  differentiation  of  F  and  G  with  respect  to  V  Yv  ,  and 
Zy  produces  the  following: 


9(F,G) 

Cy  ”  3(X  ,Y  ,Z  ) 

v  V  V 


1  /  T  0  -R  \ 

^  (  K 

\0  T  -S/ 


If  the  fundamental  angles  of  ^GVp  are  roll  (w),  pitch  ($),  and  yaw  (K), 
then  a  slightly  more  complicated  exercise  in  partial  differentiation 
produces  the  following: 


IV, .S l 

3(w,*,K) 


3F 

3F 

3F 

7w 

"5k 

3G 

3G 

3G 

7* 

■5k 

where 


S/T 

3G 

IK 

= 

Sec  $ 

_  3F  ^ 

(V 

T 

JSo 

3G 

J'V 

■SkJ 

Sec  ♦ 

_  3G 

(V 

JSo 

3F  1 

_T^v' 

■  .1  1 

T 

w\ 

3F  3F  3F 

Iw  “  (YP_YV)TZV  “  ^  ZP_ZV)'3YV 


3G  v3G  3G 

-  ■  Y  -Y  C  -  (z  -Z  W- 
3w  P  V  3Zy  K  ?  v'lY 


Define  the  (3x3)  diagonal  matrix  where  the  diagonal  elements  are 
(L+l)  dimensional  row  vectors  defined  as  (1,  t^,  t2^  . ..,  tLJK)  =  Tj. 
Then,  from  equations  (A3)  and  (A4)  ,  the  partial  derivatives  of  F  and  G 
with  respect  to  the  position  parameters  and  with  respect  to  the  attitude 
parameters  are  CV.EM  and  Ca.EN  ,  respectively.  Let  AA  be  the  3(M+1) 
vector  of  position  corrections  and  let  AB  be  the  3( N+l )  vector  of  attitude 
corrections.  The  first  (M+l)  corrections  of  AA  pertain  to  Xy,  the  next 
(ttfl)  corrections  pertain  to  Y v,  and  finally  the  last  (M+l)  corrections 
pertain  to  Zy.  A  similar  definition  holds  for  AB.  The  "linearized"  part 
of  F  and  G  that  pertains  to  the  exterior  orientation  parameter  corrections 
is  the  following: 


WA  +  C„VB 


or 


O-  c°e»]  p 

or 


QA  (B5) 

Q  :  2x3(»mH-2) 

A  j  3(MfW-2)xl 


Equations  (Bl)  and  (B2)  are  also  nonlinear  in  the  measured  data; 
namely,  line  (l)  and  sample  (p).  The  partial  derivatives  of  F  and  G  with 
respect  to  l  and  p  are  calculated  as  before  by  applying  the  chain  rule. 


3(F.G) 


3(F,G)  3(£,T)  3(F,G)  Vv’ 

<KC,Y)  3(  i,p)  +  S(xv,Y^,zvi  5 1 


SCXy.Yy.Zy) 


3t 

3(^»p) 


(B6) 


From  equation  (Bl) 

( 


3(F,G) 


3t 


3(F,G)  3(w,»,K)  _ 

3(w,*,K)  3t  3TTTpT 


Tan  Y  Sec  £  Tan  £ 
Sec2£ 


Sec2Y  Sec  £ 
0 


) 
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i  *  -  «.  -  W| 

| 


<  K 


and  from  equations  (A2)  and  (A6) , 


m  - 


(?L-  5p)  /  (Np_j) 
0 


Ax  :  line  spacing 


Define  tLT  -  (o,  1,  2tJR,  ....  LtJKL-1>jand  calculate  the  following 
vector  quantities: 


-AT. 


H  =  B  T 
a  N 


Then  from  equations  (B3)  and  (B4) 


3(F,G) 

VV9 


3<VVV  3t  „  „  3t 

- Jft -  •  TOTTpI  “  V,HV  *  TflT^r 


3(F,G)  3(w,4,K)  9t 

3(w,«|»,K)  *  3t  3(*,pT  LaHa  *  3(fc,p!) 


where 

ajlh)  '  <ST>  £T/(Np-1)) 

Note  that  time  does  not  vary  over  lines  within  a  swath. 

Let  D  represent  the  (2x2)  matrix  given  in  equation  (B6),  and 
let  (e,,e  )  represent  the  unknown  errors  associated  with  the  measured  line 
and  sample  values.  The  "linearized"  pair  of  observation  equations  then 


nlulQil 

The  matrix  elements  as  well  as  F  and  G  are  evaluated  at  the  estimated 
values  of  the  parameters.  Let  be  the  (2x2)  covariance  matrix 
associated  with  itt'  observational  pair,  then  the  ith  weight  matrix  is 

»i-  (ViV  )"' 


and  the  least-squares  solution  for  A  is 


Nq:  Number  of  observations 


There  are  (2NQ  -  3(N+-M+2))  degrees  of  freedom  associated  with  this 
adjustment.  The  residual  sum  of  squares  is 


and  the  covariance  matrix  associated  with  the  parameter  estimates  is  the 
following: 


APPENDIX  C.  LANDSAT  DATA  GENERATOR 


Assume  that  the  dynamic  model  described  by  equation  (Al)  is  accurate, 
then  if  P  is  known,  the  exposure  time  tp  can  be  determined  and  the 
corresponding  line  (i^)  and  sample  value  (Pp)  may  be  computed  from  tp. 
Rewrite  equation  (Al)  in  the  following  way: 


x  /l 
I  I 


M  M 

SV  GV 
J  JK 


Y  -  Y 
P  V 


Let  M  =  R  ,  then  if  (Cl)  is  expanded,  the  following  equation 
SVJ  07 JK  GSJK 

may  be  solved  for  the  exposure  time  tp  =  tj^.  For  convenience,  drop  the 
I,  J,  and  K  subscripts. 


E(t)  -  r  (X  -X  )  +  r  (Y  -Y  )  +  r  (Z  -Z  )  -  0 
21  p  v  22  p  v  23  p  v 


Note  that  (rjj,  *22*  r23^  ls  t'ie  second  row  of  ®GS  and  is  therefore  the 
scanner  y-axis  expressed  in  (G). 

Suppose  that  an  initial  estimate,  t  ,  of  the  exposure  time  exists, 

P0 

then  t  given  below  is  a  better  estimate: 

P1 

tp  «  t  -  E(t  )/E(t  ) 

P1  0  P0  0 


The  expression  E  refers  to  the  time  derivative  of  E.  The  process  can  be 
Iterated  until  the  correction  is  insignificant. 

The  time  derivative  is  calculated  in  the  following  way.  Rewrite 
equation  (2)  in  vector  format. 


E(  t)  *  R2  •  L 
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R2T  =  (r  ,  r  ,  r  ) 
21  22  23 


lT  -  Eyv  <W  (zp-zv>] 


The  time  derivative  of  E  then  is 


E(t)  —  R2  .  L  +  Rp  .  L 


where 


1  *  'V  V 

From  equation  (A4), 


<V  V  ZV}  =  TMT  ^ 


(Xy,  Yy,  Zy) 


•  T  T 
tma 


where 


=  (0,  1,  2tp,  Mtp 


The  time  derivative  of  the  vector  R2  can  be  computed  in  a  similar 
manner.  From  equation  (A5)  and  equation  (Bl), 


m 

m  ' 

21 

31 

m 

m 

22 

32 

m 

m 

23 

33 

Cos 
Sin  5, 


**  ’  •,  *.'  • -i,  V  v  V’ 


and 


m  m 
21  31’ 


'Cos  £T 


m  ra 
21  31^ 


-Sin  £, 


2  =  I  m  m 
1  22  32 


Sin  £t 


ra  ra 
22  32 


Cos  £. 


:) 


m  m 
23  33 1 


m  m 
23  33 , 


?p  -  (5l  "V/6t 


5P  =  ?F  +  ?P  (tP  _t 1K^ 


.th 


T  :  Start  time  of  K  swath. 
I K. 


The  time  derivatives  of  the  6  are  obtained  by  a  straightforward 
differentiation  of  the  elements  with  respect  to  time. 
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m 

31 
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K 
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22  I 

=  I  -m 
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m 

32 
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ra 

33 
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K 

-m 

13 
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m 
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33 

31 
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v  32 
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From  equation  (A3), 


(w,*,K)  =  TnT  BT 
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Recall  from  appendix  A  that  six  sample  values  are  collected 
simultaneously  over  a  swath  and  that  Xj  Identifies  the  specific  line 
within  a  swath.  From  equation  (Cl), 


x  /l 
I  I 


.  L/d  -  a 


where 


Rj  is  the  first  row  of  R^g 


Then  Xj  *  fa/Vl-a^,  where  if  f  Is  given  In  micrometers 


-3  <  xx/71  <  3 


The  line  number  within  the  swath  Is  the  nearest  integer  associated  with 
Xj/71.  The  sample  number  is  the  nearest  integer  associated  with  the 
following  expression. 


(t  -  t  -  (K-l) AT)(NP-1) 

- St - +1 


The  swath  number  K  must  satisfy  the  following  relation: 


K-l  <  (tp-  t^/AT  <  K 


Note  that  the  initial  estimate  ^Pq  must  be  within  the  correct  swath 
time  limits.  Since  the  sweep  is  normal  to  the  vehicle  path,  an  estimate 
of  the  exposure  time  can  be  obtained  by  solving  the  following  equation  for 
time: 


% *»  J  »  #  ►  m*  , ■*  m  m**  ■  •  4  •  •  ■  m  •  ***»  »•  ■.*»*.**  *  »  *  •  -  •  *  i  ‘ 


*  *  *  •  •  •'  *  ’  •  *  »  '  •  *  k*  » *  • 


G(T)  =  VVV  +  VW  +  VW  =  ° 

This  equation  can  be  solved  for  tpp  in  a  manner  similar  to  the  solution  of 
E.  In  this  case,  the  midtime  of  the  imaging  event,  say  tQQ ,  can  be  used 
as  a  first  estimate.  An  improved  estimate  is 


t01  t00  G^Too^G^To(P 


where 


G(T)  -  yXp-Xy)  +  ?v(lfp-Yv)  +  Zv(Zp-Zv)  -  «/  +YV2  +ZV2) 


where 


••  *•  ••  ••mm 

(XV,YV,ZV)  =  X  A 

X  T  -  /o,  0,  2,  6tp,  I2tp2, 


M(M-I)tpM~2 
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APPENDIX  D.  TWO-POINT  ERROR  PROPAGATION 


S- 


Consider  the  intersection  equations  expressed  by  equation  (Al).  The 
objective  here  is  to  propagate  position  and  attitude  model  errors  into  a  6 
x  6  convariance  matrix  for  two  points  Pj  and  P2*  The  right  side  of 
equation  (Al)  must  be  differentiated  with  respect  to  the  3(M+1)  position 
parameters  AA  and  with  respect  to  the  3(Nfl)  attitude  parameters  AB. 
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From  equation  (Al), 


3(X  ,Y  ,Z  ) 
PPP 

- I 

3(X  ,Y  ,Z  )  3 

VVV 


3(X  ,Y  ,Z  )  /  P 

L-  r  »  M  /  B 

3XP  GV  P 


r~  ■)  1/2 

Since  dp  -  (Xp-Xy)2  +  (Yp-Yy) 2  +  (Zp-Zy)2  I  ,  then 
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The  following  expression  represents  the  general  (i,j)c  term  of  the 
required  transformation: 


d  ^  .  <P 

y  ■  w,  K  when  j  =  1,  2,  or  3 

The  derivatives  of  the  vectors  are  determined  from  the  basic  rotation 
matrix  description  and  are  given  as  follows: 
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From  the  discussion  in  appendix  B  concerning  equation  (A4), 


3B 


and  finally, 


3(X  ,Y  ,Z  ) 

P  P  P 

-  =  WE 

3B  N 


(D2) 


3ci 

where  wtj  =  dp  -3^ 


Let  ApT  «  (AXp.AYp.AZp),  then 


4P '  (q*,Qb  ’  (“)  (D3> 

Qa:  3  x  3(MH)  matrix  of  coefficients  defined  by  equation  (Dl) 

Q  :  3  x  3(N+1)  matrix  of  coefficients  defined  by  equation  (D2) 
AA:  3(MH)  vector  of  model  position  errors 
AB:  3(NH)  vector  of  model  attitude  errors 


Fro*  equation  (D3),  the  linear  relation  between  exterior  orientation 
errors  and  errors  in  the  two  computed  ground  points  Pj  and  P2  is 
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and  from  equation  (B7),  the  (6  x  6)  covariance  matrix  for  and  P2  is 


Qa  Qb 
A1  B1 


qa 
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Qp  QD 

B1  B2 


Equation  (D4)  represents  errors  in  Pj  and  P2  due  solely  to  exterior 
orientation  errors.  Assume  that  the  ground  coordinate  frame  (G)  pertains 
to  a  DTM,  then  the  (4  x  4)  covariance  matrix  aJ  p  /  represents  the 

1  2 

horizontal  error  of  two  points  intersected  in  (G).  The  covariance 
matrix  o'  p  '  is  obtained  from  op  p  by  deleting  the  3rd  and  6  th  rows 
1  2  „  12, 
and  also  by  deleting  the  3rd  and  6th  columns. 
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APPENDIX  E.  GRAY  SHADE  RELIEF  IMAGE 


A  gray  shade  relief  presentation  is  a  mapping  that  converts  regularly 
spaced  elevation  data  into  a  terrain  surface  that  is  recognizable  as 
such.  The  example  derived  in  this  Research  Note  used  the  Lommel-Seeliger 
model  of  a  Lambertian  surface  to  produce  reflectance  values  for  each 
elevation  point. ^  A  pseudoparallax  is  introduced  so  that  a  stereomate  can 
be  produced  for  stereo  viewing. 


Consider  figure  El  where  S  is  the  direction  to  the  fictitious 
sun,  V  is  the  direction  to  the  viewer,  and  Np  is  normal  to  the  slope  at 
any  point  P  in  the  terrain  model: 


h 


Figure  El.  Reflectance  Geoaetry. 


*B.  Horn,  "Hill-Shading  and  the  Reflectance  Map,"  Proceedings:  Image 
Understanding  Workshop,  Sponsored  by:  Information  Processing  Techniques 
Office,  Defense  Advanced  Research  Projects  Agency,  April  1979,  pp.  79  to  120 


Directions  to  the  fictitious  sun  and  to  the  viewer  are  defined  for 
convenience.  For  example,  the  azimuth  and  elevation  to  the  sun  in  the 
horizontal  plane  can  be  selected  to  highlight  a  series  of  mountain 
ridges.  The  observer  can  be  placed  so  that  the  resulting  pseudoimage 
provides  a  useful  perspective.  In  order  to  prevent  line-of-sight 
ambiguities,  V,  in  the  examples  given  in  this  Research  Note,  is 
perpendicular  to  the  horizontal  plane.  The  resulting  pseudoimage  is  an 
orthographic  representation  of  the  terrain. 

The  Lommel-Seeliger  reflectance  model  is  defined  as 


The  slope  vector  Np  is  estimated  from  the  digital  terrain  matrix.  For 
example,  if  p  refers  to  the  (i,j)fc^  location  of  the  DTM,  then 
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NP  3  (~a,  -b,  1) 


(hm.j  -hi-i,3>/24Y 
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(AX, Ay):  horizontal  spacing  of  X  and  Y  in  the  DTM. 


"i.j+l  "hi,j-l)/2AX 


Valid  reflectance  values  ((XF.p<_l)  are  quantized  to  256  gray  shades  for  TV 
viewing . 


Assume  that  a  gray  shade  relief  image  of  the  terrain  is  constructed, 
one  line  at  a  time,  according  to  the  method  given  above.  The  stereomate 
is  constructed  by  shifting  each  of  its  pixel  locations  along  the  line 
according  to  computed  parallax  values.  Parallax  is  estimated  by  regarding 


the  two  pseudoimages  as  vertical  frame  images.  Consider  the  point  on 
the  line,  then  if  there  is  a  change  in  elevation  between  the  j  and 
(j+l)St  points,  the  parallax  is 


A pv  (B/H)S  (hj-hj+1) 
c 

B/H  >  0:  Base-Height 

5  =  D  /D  :  Scale 
c  c  G 

D  :  Length  of  TV  line 
c 

D  :  5116 
G 

6  =  Ax:  Horizontal  ground  spacing 

Data  is  processed  one  line  at  a  time.  The  objective  is  to  compute  a 
vector  of  coordinate  values  for  the  stereomate  line  that  corresponds  to  a 
regularly  spaced  gray  shade  relief  line  of  data.  The  computed  coordinates 
will  reflect  the  induced  parallax.  This  operation  will  produce  a  set  of 
irregularly  spaced  coordinates  that  correspond  to  the  gray  shades  of  the 
first  line.  Regularly  spaced  gray  shades  of  the  stereomate  line  are 
determined  by  interpolation. 

Let  Q  *  Sc  B/H,  then  the  Jth  coordinate  of  the  first  line  is  X,  .=  j  - 
Q(h..-hj.)  where  X.j  =1.  If  the  first  coordinate  of  the  second  line  is 
shitted  by  -QCh-.-njj),  then  X2^  *  j  -  Q(H2J-hjj)  and  in  general  Xj.  =  j  - 
Q(hij-hn).  Reflectance  values J pertaining  to  regularly  spaced  integer 
coordinates  must  be  interpolated  for  in  the  irregularly  spaced  X- 
coordinates.  The  reflectance  values  used  in  the  interpolation  are  those 
derived  for  the  first  image.  Since,  by  construction,  there  is  no  y- 
parallax  in  the  derived  model,  the  two  pseudoimages  may  be  viewed  quite 
readily  in  stereo,  for  example  as  an  anaglyph. 


APPENDIX  F.  ANAGLYPH  REPRESENTATION 


An  anaglyph  is  a  special  form  of  a  stereogram,  which  in  turn  is  a 
three-dimensional  image  of  the  terrain  derived  from  overlapping 
photographs  using  principles  of  stereoscopy.  Consider  the  case  where 
vertical  and  overlapping  photographs  are  placed  in  a  binocular  device  and 
in  such  a  way  that  the  observer  can  view  the  images  concurrently,  one  with 
each  eye.  If  the  images  are  in  the  same  plane  and  if  the  ocular  base  of 
the  viewing  device  is  parallel  to  the  photographic  base,  then  the  images 
will  fuse  and  produce  a  miniature  model  of  the  terrain.  In  the  example 
described  here,  there  is  no  y-parallax,  which  is  to  say  corresponding 
lines  of  imagery  are  parallel  to  one  another  and  to  the  ocular  and 
photographic  bases.  This  is  exactly  the  situation  between  a  gray  shade 
relief  image  and  its  stereomate  (see  appendix  E) . 

An  anaglyph  can  be  formed  by  projecting  overlapping  images  in 
complementary  colors  onto  a  common  plane .  The  anaglyph  is  viewed  by  an 
observer  using  glasses,  where  one  lens  is  the  first  complementary  color 
(such  as  red)  and  the  second  lens  is  the  second  complementary  color  (such 
as  blue).  A  color  TV  display  device  can  be  used  to  construct  anaglyphs 
from  overlapping  digital  images  or  from  a  gray  shade  relief  image  and  its 
stereomate.  In  the  first  case,  the  stereo  pair  may  need  to  be  reformatted 
to  remove  y-parallax.  In  either  case,  the  images  are  displayed 
simultaneously  and  in  complementary  colors.  The  anaglyph  is  viewed  by  the 
observer  using  the  appropriate  glasses. 

The  measuring  techniques  used  in  analog  operations  can  be  applied  in 
the  digital  realm.  For  example,  consider  the  situation  where  overlapping 
images  are  viewed  as  an  anaglyph.  Digital  image  marks  defining  the 
line  of  sight  for  each  eye  can  be  superimposed  on  the  anaglyph  in 
complementary  colors.  If  the  marks  are  placed  over  corresponding  image 
points,  then  the  two  will  fuse  and  appear  to  lie  on  the  image  terrain. 
Since  the  location  of  the  marks  on  the  TV  displays  are  uniquely  related  to 
the  data  base,  the  described  operation  is  similar  to  the  operation  of  a 
stereocomparator.  In  much  of  the  digital  image  analysis  at  the  Computer 
Sciences  Laboratory,  ETL,  corresponding  points  are  defined  to  be  at 
regular  line  and  sample  intervals  on  the  left  image.  Thus,  if  there  is 
little  or  no  y-parallax,  the  cursor  on  the  right  image  can  be  moved  along 
a  TV  line  to  affect  correspondence.  An  application  of  these  principles 
was  tested  and  described  more  completely  by  ETL's  Norvelle.^ 


2 

F.  Norvelle,  Interactive  Digital  Correlation  Techniques  for  Automatic 
Compilation  of  Elevation  Data,  U.S.  Army  Engineer  Topographic  Laboratories, 
Fort  Belvoir ,  VA,  October  1981,  ETL-0272,  AD-A109  145. 
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